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Abstract. A mean-field approach (MFA) is proposed for the analysis of orien- 
tational order in a two-dimensional system of stochastic self-propelled particles 
interacting by local velocity alignment mechanism. The treatment is applied to 
the cases of ferromagnetic (F) and liquid-crystal (LC) alignment. In both cases, 
MFA yields a second order phase transition for a critical noise strength and a scal- 
ing exponent of 1/2 for the respective order parameters. We find that the critical 
noise amplitude r) c at which orientational order emerges in the LC case is smaller 
than in the F-alignment case, i.e. rj^ c < T}q. A comparison with simulations of 
individual-based models with F- resp. LC-alignment shows that the predictions 
about the critical behavior and the qualitative relation between the respective 
critical noise amplitudes are correct. 

1 Introduction 

Self-propelled particles (SPPs) with local interactions can organize into large-scale patterns. 
Flocks of birds [1], swarms of bacteria [213] . sperm cells [4], mixtures of microtubules and 
motors [5], are examples for such a behavior. Beyond the complexity of each particular system, 
we observe that there are some few common features which cause the emergence of long-range 
order in these systems: the active motion of the particles and a velocity alignment mechanism. 

The Vicsek-model [6] is considered the simplest model of SPPs which exhibits collective 
motion. In this model, point-like particles moving with a velocity vector of constant magnitude 
interact by aligning their velocity direction to the local average velocity. One can think of this 
model as a model of moving spins, in which the velocity of the particles is given by the spin- 
vector. Going further in this analogy with spin systems we denote this alignment mechanism 
as ferromagnetic (F-alignment). The temperature associated with spin-systems enters in the 
Vicsek-model as noise in the alignment mechanism. It was shown that two-dimensional SPPs 
with F-alignment and additive noise exhibit a second-order phase transition which leads to low 
enough noise to long-range orientational order [61718] . For different choices of system parameters 
and different choices of the noise term (non-additive noise term), simulations show, however, 
first order phase transitions to orientational order [5], Interestingly, in equilibrium systems of 
non-moving spins with continuum symmetry this transition cannot occur [10] , 

F-alignment is one possible alignment mechanism, but clearly not the only one. If a system of 
self-propelled rods interacts simply by volume exclusion as described in [llj , particles may end 
up moving in the same direction as well as in opposite directions. A biological realization of such 
a system are myxobacteria, which in the early stages of their life-cycle organize their motion by 
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simply pushing each other [3]. A similar effect without active motion occurs in liquid-crystals 
at high density where particles get locally aligned [12]. In analogy to these systems we name 
this mechanism hence liquid-crystal alignment (LC- alignment). In a system of SPPs with LC- 
alignment particles align their velocity to the local average director. In simulations of a model 
of SPPs with LC-alignment it was found that at high density these moving "liquid-crystal" 
spins exhibit a second-order phase transition leading to long-range orientational order for low 
noise [13] . Notice that the orientational order observed in SPPs with LC-alignment refers to the 
emergence of a global director in the system, while for F-alignment orientational order refers 
to the appearance of a global direction of motion. 

Toner and Tu were the first to look for a macroscopic description of SPPs with F-alignment. 
Their approach was a phenomenological hydrodynamical description based on symmetry argu- 
ments for which they derived general macroscopic equations for a large class of individual-based 
models of SPPs with F-alignment [14I15|16] , In this approach many of the parameters in the 
model are difficult to derive from the microscopic dynamics. Recently Gregoire et al. proposed 
an alternative approach based on the Boltzmann equation and could explain in a systematic 
way the functional form proposed by Toner and Tu [17j by use of an ad-hoc collision term. 

The case of LC-alignment has been much less explored. Recently, Ramaswamy et al. pro- 
posed a phenomenological hydrodynamical description for driven, but non-persistent, particles 
with LC-alignment [18] , The approach is comparable to the Toner and Tu description for F- 
alignment. One of the striking results of this approach is the giant number fluctuations of 
particles in the ordered state, which has been confirmed in simulations by Chate et al. [19] . 
More recent simulations have suggested that these fluctuations are linked with intrinsic phase 
separation of SPPs into regions of high and low density [20] . 

Here, a mean-field type description for SPPs with F- and LC- alignment is derived. Nu- 
merical evidence provided by individual-based simulations indicates that SPPs with both F 
and LC-alignment (and additive noise) can exhibit a continuous kinetic phase transition in 
two dimensions. The derived mean-field equations allow us to study ferromagnetic as well as 
liquid-crystal interactions among particles. Through this approach the phase transition to orien- 
tational order observed in individual-based simulations at high density for F and LC-alignment 
is correctly captured. In addition, we show that the critical noise amplitude r\ c is such that 
Vc C < m the mean-field description as well as in the individual-based simulations in both 
investigated cases. 

The paper is organized as follows. In section 2 we introduce an individual-based model 
of SPPs which interact by either F- or LC-alignment, and give the corresponding definitions 
of the order parameters. In section 3 we present a macroscopic description of the individual- 
based model introduced in Section 2. We simplified the corresponding macroscopic model by 
considering a mean-field ansatz, and perform a linear stability analysis of the equations. In 
section 4 we compare the mean-field description with high density simulations in the limit of 
very fast angular relaxation. We discuss the limitations of the mean-field approach in section 
5. 



2 Individual-based model 
2.1 Equations of motion 

We consider point-like particles moving at constant speed in a two dimensional space and 
assume an over-damped situation such that the state of particle i at time t is given by its 
position Xi and its direction of motion 6i. The evolution of these quantities follow: 



Xj = v o v{0i) 



(1) 




7— (x h 9i)+fji(t) 



(2) 



where 7 is a relaxation constant, and U the interaction potential between particles, and hence 
Qi) defines the velocity alignment mechanism. Moreover, v represents the active velocity 
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Fig. 1. Temporal evolution of the velocity direction distribution (angular distribution) in simulations 
with very fast angular relaxation, (a) corresponds to F-alignment, while (b) to LC-alignment. Number 
of particles TV = 100, radius of interaction e = 2, linear system size L = 42.4, and noise amplitude 
77 = 0.25. 



of the particles, v(0i) is defined as v(#i) = (cos(#i), sin(#i)). The noise applied to the direction of 
motion, fji(t), obeys the following statistics: (rji(t)) = and {fii{t)fij(t')) = rS(t — t')Sij, where 
8[t — t') is a Dirac-delta function, S^j is a Kronecker-delta function, and t is the "strength" of 
the noise. The evolution Eqs. ((T|) and @ are expressed in terms of first derivatives. In this way, 
vo in Eq. (fT|) can be considered as an active force divided by a translational friction coefficient, 
and 7 in Eq. as the inverse of a rotational friction coefficient. 

In analogy to spin systems, the ferromagnetic velocity alignment mechanism is given by a 
potential defined as: 

U F (x l ,6 i ) = - <x»(0i-0j) (3) 

l x i — X J I ^ e 

where e is the radius of interaction of the particles. For the liquid-crystal alignment mechanism, 
we choose the potential introduced by Lebwohl and Lasher to study liquid crystal interactions 
on a lattice [3T] which reads: 

U LC {xiA) = - cos 2 ^-^) (4) 

One can add a coupling strength coefficient to the expression ([3]) and ((4]). We assume that 
the coupling strength is absorbed in 7 in Eq. Notice that the potential given by Eq. ([3]) 
exhibits one minimum, while Eq. (fj[|) has two minima, which correspond to particles pointing 
in the same direction and particles pointing in opposite directions. 

In the limiting case of very fast angular relaxation we obtain from Eqs. (p} and |2]) the 
updating rules: 



n +At - arg 



v Q v (91) At (5) 

Y f(v(0'.),v(0i)) )+vi (6) 

where arg (b) indicates the angle of a vector b in polar coordinates, and r\\ is obtained from a 
distribution p(r]l) defined as p(rjf) = l/i] when 77* belongs to the interval [— f , §] , and p(r)\) = 
otherwise. In consequence, (r)\) = and (f?f?7j ) = (rj 2 /12)Sij8t,t' ■ Given two vectors a and b, 
the function f(a, b) is defined as follows. For F-alignment, f(a, b) = a. For LC-alignment, f 
takes the form: 

a if a.b > 

-a if a.b < ^ > 
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as described in [13]. Notice that F-alignment implies that the sum in Eq. ((6]) becomes a simple 
weighted local velocity average. In consequence, under F-alignment Eqs. ((5]) and (J6j> define 
the original Vicsek model [6]. Particles interact by calculating the local average direction of 
motion, and if the noise strength is low enough move roughly in that direction. In contrast, 
LC-alignment implies a local average of mapped velocities that leads particles to calculate the 
local average director (and not the average direction of motion as for F-alignment). This process 
defines locally two possible directions of motion, and particles choose from these two options 
the one that is closer to their present direction of motion. Notice that throughout the text F- 
alignment refers to the Vicsek model, and this can be either in its original discrete-time form, 
or in its generalized continuum time form given by Eq. fl]), (0, and ([3]). 



2.2 Order parameters 



If particles interact through the F-alignment mechanism, and assuming low noise amplitude, 
they get locally aligned, and locally point in a similar direction. The question is whether such 
local alignment may lead to a global orientational order in which a macroscopic fraction of the 
particles in the system points in a similar direction. The order parameter that quantifies this 
phenomenon is the modulus of the normalized total momentum (analogous to the magnetization 
in the XY-model|10j) that we express as: 



S* = 



1 N 



i=0 



(8) 



where N stands for the total number of particles in the system. S F takes the value 1 when all 
particle move in the same direction. On the other hand, S F is equal to in the disordered case in 
which particles point in any direction with equal probability. This can be also observed through 
the velocity direction distribution, that in two dimensions becomes an angular distribution C{6). 
For high values of the noise, C(9) is flat. When the noise is decreased below a critical noise 
r\ c an instability arises in the system (characterized by a single peak) indicating the onset of 
orientational order as shown in Fig.QJa). 

On the other hand, if for example, half of the particles move in one direction, and the 
other half in the opposite direction, S F is also 0. Clearly, S cannot distinguish such a state 
and the completely disordered state. However, LC-alignment may induce such a kind of local 
arrangement of particle velocities, and lead to a global orientational order state in which there 
are two opposite main directions of motion in the system. To study such orientation ordering, 
one uses the order matrix Q of liquid crystals [12] ■ For two dimensions one takes the largest 
eigenvalue S LC of Q and obtain the following scalar orientational order parameter: 



S LC = 



\ 




Vxi'VyiVxjVyj 



(9) 



where 



and 



are defined as 



cos(6*i) and 



'yl 



sin(#j). The orientational order 



parameter S LC takes the value 1 when all particles are aligned along the same director, and the 
value i in the disordered phase where particles move with equal probability in any direction. 
Again this can be observed through the velocity direction distribution C(9). In this case, for 
low values of the noise amplitude, as shown in Fig. [lib), an instability arises in the system with 
the characteristic of having two peaks separated by 2tt. 
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3 Mean-field approach 

3.1 Derivation of the mean-field approach 

A system of SPPs may alternatively be described through a density field ip(x,v(9),t) = 
ip(x, 9, t) in such a way that the particle density at a point x is given by 

p(x,t) = / ip(x,6,t)d0 (10) 
Jo 

while the velocity direction distribution (or angular distribution) can be expressed as: 

C(0,t) = [ ip(x,6,t)dx (11) 
Jn 

We recall that in the individual-based model the kinetic energy is conserved, while the 
momentum is not. For F-alignment, the system tends to increase the total momentum, while 
for LC-alignment the tendency is to decrease it. The continuum approach has to reflect that 
particles can re-orient their velocity direction but always move at constant speed. On the other 
hand, the number of particles has to be conserved. Under these assumptions the following 
evolution equation for ip(x,9,t) is obtained: 

d t 4> = Dgdeeip - d 9 [F e tp] - V M (12) 

where Fgip and F x ijj are deterministic fluxes which are associated to the local alignment mech- 
anism and active migration, respectively, and Dg refers to the diffusion in the direction of 
motion. 

Let us derive the specific expressions for Dg, Fg and F x . Dg depends on the square of 
the noise amplitude. For example, in the individual-based simulations rji(t) has been taken, as 
mentioned above, from a homogeneous distribution of width rj and centered around 0, and in 
consequence Dg is given Dg = rj 2 At/ '24, where At is the temporal time step. Fg contains the 
interaction of a particle located at x and pointing in direction 9 with all neighboring particles 
which are at a distance less than e from x, and so takes the form: 

Fg = - 7 / dx' /%' ^fo g ' x/ '% (a; ',fl',t) (13) 

Jfl(x) JO 0U 

where U(x, 9, x', 9') represents the pair potential between a particle located at x and pointing 
in direction 9, and another at x' and pointing in direction 9'. R(x) denotes the interaction 
neighborhood around x. If x' is inside R(x), then U(x,6,x' ,9') = U(9,9'). Finally, Fg can 
be thought as the "torque" felt by a particle located at x and pointing in direction 9. The 
expression for is straightforward and is directly related to the velocity of particle at x and 
pointing in direction 9, 

F x = v v(9) (14) 



3.2 Angular distribution 



Integrating both sides of Eq. (|12[) over the space Q we obtain an evolution equation for C(9, t) 
which still depends on ip{x, 9, t). In the following we assume a homogeneous spatial distribution 
of particles ip(x,9,t) = C(9,t)po/N, where po is defined as p — N/L 2 , being L the linear size 
of the system. With these assumptions the equation for the temporal evolution of C(9, t) reads: 



dC{9,t) 
dt 



= D, 



d 2 C{9,t) 
89 2 



+ l-prdg 



o 



C(9,t) 



(15) 
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(a) 




(b) 



100 




100 



Fig. 2. Temporal evolution of C(8,t). (a) F-alignment, numerical integration of Eq. (|16[) with Dg = 
0.28. (b) LC-alignment, numerical integration of Eq. ([20|) with Do = 0.014. For both (a) and (b), 
C* — 0.3183, At = 0.001 and AO — 0.16. The initial condition is a random perturbation around C*. 
Notice that for F-alignment a single peak emerges, while for LC-alignment the distribution develops 
two peaks. 



3.3 Linear stability analysis for F-alignment 

For both F- and LC-alignment the homogeneous angular distribution is a steady state of Eq. 
(fl~5)) . We determine the onset of the ordered state by studying the linear stability of the disor- 
dered state. First let us look at the F-alignment. By dividing both sides of Eq. (fl5|) by "fire 2 /L 2 , 
and redefining time as r = ("/ire 2 / L 2 )t, and D' g — Dg /["fire 2 / L 2 ] one obtains: 



dC (0,t) 
JTt 



D' e dggC{e,t)+dg 



<W'sm(6-0')C(9',t) >C(6,t) 



Now, consider a weak perturbation of the homogeneous pattern: 

C{6,t) = C* +C e m9 e XT 



(16) 



(17) 
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Notice that e are eigenfunctions of the operators emerging from the linearization of Eq. (fl6|) 
about the homogeneous steady state. By substituting into Eq. (|16j) and keeping terms linear in 
Co we obtain the following expression for the eigenvalues: 

Re(X) = -D' e n 2 + nC*6 nA (18) 

This means that the only mode which can become unstable is n = 1. The condition for the 
instability of the homogeneous state takes the form: 

P0 > ^§ (19) 

where po = N/L 2 . For a given noise amplitude, expressed by Dg, there is a critical particle 
density above which the homogeneous solution is no longer stable. Fig. QJa) shows that in 
individual-based simulations indeed a single peak emerges in the system for low density. Fig. 
[2]^a) confirms that such qualitative behavior is recovered by numerical integration of Eq. (1161) . 



3.4 Linear stability analysis for LC-alignment 

Applying analogous procedure for LC-alignment yields: 
dC(9,t) 



Or 



D e deeC(6,t)+d e 



dd'2 cos (9 - 6') sin (0 - 9') C (9', t) } C {9, t) 



(20) 



Again the weakly perturbed homogeneous ansatz given by Eq. (TIT)) is considered. As before e m6 
are eigenfunctions of the linearized operators. Substituting Eq. (fT7|) into Eq. (|20|) and keeping 
terms linear in Co the following expression for the eigenvalues is obtained: 

Re(X) = ~D' g n 2 + 2TrC*S n , 2 (21) 

As for the F-alignment, there is only one mode which could become unstable, but this time it is 
n = 2. This mode is the only one that exhibits two peaks separated by n, which corresponds to 
two population of particles migrating in opposite direction. In this case, the instability condition 
of the homogeneous states takes simply the form: 

Po > ~~2 ( 22 ) 

Again, this inequality defines a critical density for a given noise amplitude above which the 
homogeneous solution is no longer stable. Fig.QJb) shows the emergence of these two peaks for 
LC-alignment in individual-based simulations. Numerical integration of Eq. (|2U|) . see Fig. fJKb), 
confirms that this behavior is recovered qualitatively by the mean-field description. 

Eqs. (fl9|l and (|22|) indicate that the instability of the homogeneous state is given by po, Dg, 
and e, the range of interaction. The critical density is inversely proportional to e 2 , hence when 
e goes to infinity the critical density goes to 0. The interpretation of this is straightforward, 
e — > oo indicates that particles have infinity "visibility", i.e., each particle can sense all other 
particles in the system. In this way, the collective behavior has to emerge independent of particle 
density. The other limiting case is represented by e — > 0. In this case particles do not interact 
and in consequence no organized motion is possible. 

From these findings a phase diagram is derived that shows where the system exhibits velocity 
orientational order (see Fig.[3|). 



3.5 Spatially inhomogeneous steady states 



Through the linear stability analysis it has been found for which conditions the homogeneous 
distribution (disordered state) becomes unstable. To study the nonlinear behavior of these 
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Fig. 3. Phase diagram derived from the continuum approach. The unstable region corresponds to the 
velocity orientational order, while stable means no orientational order can be observed. 




3.14 Q 6.28 3.14 Q 6.28 



Fig. 4. Convergence towards the non-trivial stable steady state, (a) F-alignment, numerical integration 
of Eq. |16J) with D e = 0.3472. (b) LC-alignment, numerical integration of Eq. (J20j) with D = 0.2813. 
For both (a) and (b), C* = 0.3183, At = 0.001 and A6 = 0.0785. The initial condition is a random 
perturbation around C*. Different curves correspond to different times. Notice that for large values of 
t curves start to overlap on top of each other. 

instabilities in more detail, Eq. (|15| can be integrated numerically. Details about the numerical 
methods are given in the Appendix. Fig. [2] has shown already the temporal evolution of C(9, t). 
The initial condition is a homogeneous state with small random perturbations: C(9 n ,t = 0) = 
C* + T](n), where 6 n denotes the discrete angular variable, C* is the constant unperturbed 
homogeneous state, which we have set to be in the unstable regime according to Eqs. (fl"9]) 
and (|22|) for the F- and LC-alignment case, respectively, and f](n) is a white noise. In Fig. 
@] the angular distribution for F- and LC-alignment is shown at different times. C(9, t) tends 
asymptotically to a non-trivial steady state, see Fig. [4j The width of the peaks in the steady 
state is the result of the competition between influence of rotational diffusion, indicated by Dg, 
and the alignment force associated with the interactions. 



3.6 Scaling of the order parameter close to the transition 

For a given density, there is a critical Dg c . Close to Dg c we expect to observe that only one 
mode dominates C(6,t). As said before, n = 1 is dominant for F-alignment and n — 2 governs 



Will be inserted by the editor 



9 



LC-alignment. The steady state distribution C s t{9) then takes the form: 

C st {9) ~C* + B 1 ^D 8c - D g cos(0 - 6 ) (23) 

for F- alignment, while for LC-alignment the expression reads: 

C st {6) ~ C* + B 2 ^D 9c - D e cos(2(0 - O )) (24) 

where B\ , B\ are constants and 8q is an arbitrary phase which depends on the initial condition. 
In both cases the maximum amplitude of C s t (9) close to the Dg a grows as ^JDg c — Dg. Inserting 
Eq. ([23| into Eq. ([8|) and using, as indicated above, Dg = rj 2 At/24, we obtain the scaling of 
the order parameter S : 

S F ~ B x ^—^ (25) 

where B is a constant. To obtain the scaling of the order parameter S , we insert (|2"3T) into 
Eq. ©: 

S LC ~ i + B 2 ^~T, (26) 
where again i?2 is a constant. Si and B 2 are constants proportional to rj c . 



4 Comparison with individual-based simulations 



Individual-based simulations have been performed in the limit case of very fast angular re- 
laxation [617] . In contrast, our mean- field description assumes that there is a finite angular 
relaxation. Can we expect the mean-field approach to describe scaling of the orientational dy- 
namics in this kind of simulations? We redefine 7 as function of the particle velocity vq and the 
particle density p. The effective resulting mean-field equation reads: 



dC_ 

dt 



(0,t)=Dg 



d 2 C(9,t) 
d0 2 



de 



,017(0,1 



00 



-c{o',t)\c{e,t) 



(27) 



where 7(1*0, p) is an effective interaction strength which absorbs the spatial dynamics. 

The scaling obtained from individual-based simulations may now be compared with the one 
predicted by the mean-field approach. We recall that Dg ~ rj 2 , where 77 is the orientational 
noise amplitude used in individual-based simulations. From this we find that Dg c — Dg has 
to be Dg a — Dg = K{rj c — 77) + 0((ri c — rj) 2 ), where if is a constant. We focus on the LC- 
alignment and replace this expression into Eqs. and (|2"6")) . We obtain that C s t(0) — C* + 
B 1 ^]~ r rj cos(2(0 - 6 )) and S LC 



77, where B\ and B 2 are constants. 
Fig. Efa) shows a comparison between the scaling predicted by the mean-field approach for 
S LC (dashed curve) and the one obtained from individual-based simulations for p = 4 in the 
limit of very fast angular relaxation (symbols). We find good agreement between the mean- 
field prediction and the simulations for the scaling of S near ?y c that suggests that individual- 
based simulations with LC-alignment at high densities exhibit a mean-field type transition. 
Notice that simulations start to deviate exactly at the point where density fluctuations become 
important (denoted by the dot-dashed vertical line in ©(a)). Let us recall that the mean-field 
approach implies the assumption of homogeneous density. Evidence also points towards a mean- 
field transition if we look at the scaling of the maximum amplitude of the angle distribution as 
function of the angular noise intensity r\ (see Fig. (0(b)). The order parameter scaling exponent 
for F-alignment (Vicsek-model) has been found to be 0.45 ± 0.07 [617) . which is also in line 
with the predictions of the presented mean-field theory. 

Finally, Fig. [S] shows that in individual-based simulations with the same parameters and 
different (namely LC- and F-) alignment mechanism, in the limit of very fast angular relaxation 
rjQ < i]£ as predicted by the mean-field theory. Note, however, that the simulations yield 
w rjQ, while the mean-field description predicts V2rj^ c = r] F . 
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Fig. 5. (a) Scaling of the scalar order parameter S LC with the noise amplitude t). (b) Scaling of the 
maxima in C s t{6) with r\. Symbols correspond to simulations with LC-alignment in the limit of very 
fast angular relaxation, p — 4 and N = 2 12 . The dashed-curve corresponds to the scaling predicted by 
the mean-field approach (see Eq. (|26[l ). The vertical dot-dashed line indicates the onset of clustering 
effects in the simulations. To the right of that line particle density p(x) can be considered constant. 
To the left of the line, p(x) becomes a function of x. For F-alignment in limit of very fast angular 
relaxation, the scaling of S F is close to the one predicted by Eq. (|25p . compare also [617] . 



5 Concluding remarks 

We have derived a mean-field theory for self-propelled particles which accounts for F- and 
LC-alignment. This approach predicts a continuous phase transition with the order parameter 
scaling with an exponent one half in both cases. In addition, it yields that the critical noise 
amplitude below which orientational order emerges is smaller for LC-alignment than for F- 
alignment, i.e., r^P < r/^. 

These findings were confirmed by individual-based simulations with F- and LC-alignment. 
In the limit of infinitely fast angular relaxation used in simulations here the mean-field theory 
provides a good qualitative description of the simulations. If simulations were performed by 
integrating Eqs. ([I]) and (J2J) with a finite angular relaxation, i. e. a finite 7, a direct corre- 
spondence between parameters used in simulations and parameters in the mean-field theory 
can be made and quantitative comparisons become possible. The presented mean-field theory 
is not an exact coarse-grained description of Eqs. |T]) and {2|). For instance, we have neglected 
the potential impact of particle-particle correlations. Furthermore, we have assumed spatial 
homogeneous density to study the emergence of orientational order. Thus, the presented ap- 
proach does not apply to situations where self-propelled particles show clustering at the onset 
of orientational order [11113) . In summary, a better understanding of the problem should imply 
the study of the interplay between local orientational order and density fluctuations. We leave 
that for future research. 
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6 Appendix - numerical integration scheme 

The numerical integration of the integro-partial differential equation (|15p requires to perform 
the numerical integration of Eq. (| 1 3 j) to then proceed to the integration of the diffusive and 
advective terms in Eq. (fT5|) . 

At each time step Fn(6,i) is calculated through a simple Newton-Cotes method. Then the 
integration of Eq. (fT5|) is performed through an operator splitting method. The diffusion is 
implemented by an explicit forward method. The integration of the active turning (Eq. (fl3|) ) 
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Fig. 6. Comparing simulations of particles with F-alignment (crosses) and LC-alignment (circles) in 
the limiting case of very fast angular relaxation. In both cases N = 2 14 and p = 2.0. Notice that the 
order parameter for F-alignment is S F while for LC-alignment is S LC (See text). The dashed horizontal 
line indicates the minimum value that S LC could take. The dashed curves correspond to the best fit 
assuming an exponent 0.5, i.e., rj c was the fitting parameter. 



contained in the advective term requires special attention. Since Fq depends explicitly on 9 and 
t neither a Lax nor an Upwind method gives a satisfactory result. We overcame this difficulty 
by implementing the following variant of the Upwind method: 

C(6 k ,tj + 1) = (1 - Wnie^t^CiOk,^) + 
9(Fn(6 k - l,t 3 ))\F n {6 k - l,tj)\C(6 k - l,t 3 ) + 

e{-F a {6 k + l,t 3 ))\F f2 (0 k + l,t 3 )\C{9 k + l,t 3 ) (28) 

where k and tj represent the discrete indices of the angular and temporal variables respec- 
tively, 0{x) denotes a Heaviside function, and FQ(9 k , tj) is defined as Fa(6 k , tj) — (At/A0)Fa(9 k ,tj) 
where A9 and At are the discretization of the space and time respectively. 
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